/***********************************************************************************
PROGRAM: 02_MW_Analysis.do
PURPOSE: Runs mortality analyses for Miller, Johnson, and Wherry
DATE: 	 October 15, 2020
NOTES: 	
***********************************************************************************/
    
global data /projects/mortality/data/
global userdata /projects/mortality/data/
global output /projects/mortality/output/


*******************************************************************************************************************
***Descriptive Statistics

use $data/panelfinalelig, clear

**Counterfactual mortality rates:

gen preyear=.
replace preyear=1 if treated14==1 & dyear==2013
replace preyear=1 if treated15==1 & dyear==2014
replace preyear=1 if treated16==1 & dyear==2015
replace preyear=1 if treated17==1 & dyear==2016

**Main sample, by event year:
mean diedpanel [pweight=pwgt] if target1==1 & age2014 > 54 & sameyear==0
mean diedpanel [pweight=pwgt] if target1==1 & age2014 > 54 & preyear==1 & sameyear==0

mean diedpanel [pweight=pwgt] if event3==1 & target1==1 & age2014 > 54 
mean diedpanel [pweight=pwgt] if event2==1 & target1==1 & age2014 > 54 
mean diedpanel [pweight=pwgt] if event1==1 & target1==1 & age2014 > 54 
mean diedpanel [pweight=pwgt] if event0==1 & target1==1 & age2014 > 54 
mean diedpanel [pweight=pwgt] if posttreated==1 & povpi > 399 & age2014 > 54 & age2014 < 65
mean diedpanel [pweight=pwgt] if posttreated==1 & age2014 > 64
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 18
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & age2014 < 62
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & blacknh==1
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & whitenh==1
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & hispanic==1
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & hispanic==0 & whitenh==0 & blacknh==0
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & msp=="1"
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & msp!="1"
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & povpi < 139
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & lths==1
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & noinsure==1
mean diedpanel [pweight=pwgt] if posttreated==1 & age2014 > 64 & target1b==1
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & female==1
mean diedpanel [pweight=pwgt] if posttreated==1 & target1==1 & age2014 > 54 & female==0
   
mean diedpanel [pweight=pwgt] if posttreated==1 & newtarget==1 & age2014>54
mean diedpanel [pweight=pwgt] if posttreated==1 & newtarget==1 & age2014 > 54 & newinc_fpl<=1.38
mean diedpanel [pweight=pwgt] if posttreated==1 & newlyelig==1 & age2014>54 

**Descriptive statistics
summ diedpanel if (treated14==1 | treated15==1 | treated16==1 | treated17==1) & sameyear==1 & target==1
summ diedpanel if target==1 & age2014 > 54 & preexp==1 & sameyear==0 
summ diedpanel if povpi > 399 & age2014 > 54 & age2014 < 65 & sameyear==0
summ diedpanel if target==1 & sameyear==0 
summ diedpanel if target==1 & sameyear==0 & age2014 > 54 
summ diedpanel if target==1 & sameyear==0 & age2014 > 54 & noinsure==1 
summ white black hispanic female noinsure medicaid povpi lths if sameyear==1 & age2014 > 54 & age2014 < 65 & target==1 

****************************************************************************************************************
***Regression Analyses

global controls black white hispanic female i.age
global econcty unemp_rate medhhinc povrate
global opioid mustaccess pdmp painlaw medmj dispmj i.triplicate##i.dyear

use $data/panelfinalelig, clear

***Simple Diff Diffs (No Controls):

**Age 55-64 (main sample)
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Baseline") append 
estat vce 

*All Ages
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, append bdec(10)

***Estimate event study regression (No Controls)
*All Ages
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs.xls, append bdec(10)

** Age 55-64 in 2014
cap noi reghdfe diedpanel eventm5 eventm4 eventm3 eventm2 eventm1 event1 event2 event3 event4 sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
*outreg2 using $output/panelregs.xls, append bdec(10)
mat v=e(V)
mat b=get(_b)'
mat list v
mat list b

**************************************************************************************
*Simple Diff Diffs (With Controls):
*Target sample
cap noi reghdfe diedpanel posttreated sameyear $controls [pweight=pwgt] if target1==1, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff_controls.xls, append 

**Target sample + Age 55-64 in 2014
cap noi reghdfe diedpanel posttreated sameyear $controls [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff_controls.xls, append 

**Target sample + Age 55-64 in 2014 + labor demand 
cap noi reghdfe diedpanel posttreated bartik_allind sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Labor Demand") append 

**Target sample + Age 55-64 in 2014 + county economic controls 
cap noi reghdfe diedpanel posttreated $econcty sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls,  ctitle("Econ Controls") append 

**Target sample + Age 55-64 in 2014 + state opioid controls 
cap noi reghdfe diedpanel posttreated $opioid sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Opioid Controls") append 

**Target sample + Age 55-64 in 2014 + Trade Controls
cap noi reghdfe diedpanel posttreated c.tradeotch2000_2014##i.dyear sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Trade Controls (Excl AK HI)") append 

**Target sample + Age 55-64 in 2014 + Demo Controls
cap noi reghdfe diedpanel posttreated $controls sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Demographic Controls") append 

**Target sample + Age 55-64 in 2014 + All Controls
cap noi reghdfe diedpanel posttreated bartik_allind $econcty $opioid $controls c.tradeotch2000_2014##i.dyear sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("All Controls (Excl AK HI)") append 

    
****************************************************************************************
*Event Studies (With Controls):
**Target sample + Age 55-64 in 2014 + demographic controls
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 $controls sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls,  ctitle("Demo Controls") append 

**Target sample + Age 55-64 in 2014 + county economic controls 
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 $econcty sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls,  ctitle("Econ Controls") append 

**Target sample + Age 55-64 in 2014 + state opioid controls 
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 $opioid sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Opioid Controls") append 

**Target sample + Age 55-64 in 2014 + Trade Controls 
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 c.tradeotch2000_2014##i.dyear sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Trade Controls (Excl AK HI)") append 

**Target sample + Age 55-64 in 2014 + All Controls
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 bartik_allind $econcty $opioid $controls c.tradeotch2000_2014##i.dyear sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("All Controls (Excl AK HI)") append 

**Labor demand
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 bartik_allind sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs.xls, ctitle("Labor Demand") append 

*
*********************************************************************************************
***Heterogeneity by Marital Status
**Target sample + Age 55-64 in 2014 + Married
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & msp=="1", absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Married") append

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & msp=="1", absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Married") append

**Target sample + Age 55-64 in 2014 + Not Married
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & msp!="1", absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("Not Married") append

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & msp!="1", absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Not Married") append

**********************************************************************************************
***Heterogeneity by Race
gen blacknh=black-hispanic
replace blacknh=0 if blacknh < 0

gen whitenh=white-hispanic
replace whitenh=0 if whitenh < 0

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & blacknh==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Black NH") append 

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & blacknh==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Black NH") append 

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & whitenh==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("White NH") append 

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & whitenh==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("White NH") append 

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & whitenh==0 & blacknh==0 & hispanic==0, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Other NH") append 

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & whitenh==0 & blacknh==0 & hispanic==0, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Other NH") append 

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & hispanic==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Hispanic") append 

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & hispanic==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Hispanic") append 

**********************************************************************************************
***Heterogeneity by Gender
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & female==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Women") append

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & female==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Women") append

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & female==0, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Men") append

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & female==0, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Men") append

*********************************************************************************************
***Other Subgroups
**Target sample + Age 55-64 + LT HS only in 2014
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & lths==1 & age2014 > 54, absorb(st dyear year) cluster(st)
mean died [pweight=pwgt] if e(sample)
outreg2 using $output/diffdiff.xls, ctitle("LT High School") append

**Target sample + Age 55-64 + pov rate < 138 only in 2014
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & povpi < 139 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, ctitle("LT 138")  append
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & povpi < 139 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("LT 138")  append

**Target sample using SHADAC family income defn + Age 55-64 in 2014
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if newtarget==1 & age2014 > 54, ///
absorb(st dyear year) cluster(st)
outreg2 using $output/shadac.xls, ctitle("Target") append 
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 ///
sameyear [pweight=pwgt] if newtarget==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/shadac.xls, ctitle("Target")  append
    
**Target sample + Never Mcare Elig (55-61)
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & age2014 < 62, absorb(st dyear year) cluster(st)
mean died [pweight=pwgt] if e(sample)
    
**Target sample + Never Mcare Elig (55-61)
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & age2014 < 62, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, append 
mean diedpanel [pweight=pwgt] if target1==1 & age2014 > 54 & age2014 < 62 & posttreated==1

**Target sample + Age 55-64 in 2014 + not insured
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & noinsure==1, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, append 
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 & noinsure==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs.xls, append bdec(10)

**3-year age bins
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 73 & age2014 < 77, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("74-76") append 
mean diedpanel [pweight=pwgt] if target2==1 & age2014 > 73 & age2014 < 77 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 70 & age2014 < 74, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("72-73") append 
mean diedpanel [pweight=pwgt] if target2==1 & age2014 > 70 & age2014 < 74 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 67 & age2014 < 71, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("68-70") append 
mean diedpanel [pweight=pwgt] if target2==1 & age2014 > 67 & age2014 < 71 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 64 & age2014 < 68, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("45-54") append 
mean diedpanel [pweight=pwgt] if target2==1 & age2014 > 64 & age2014 < 68 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 61 & age2014 < 65, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("62-64") append 
mean diedpanel [pweight=pwgt] if target2==1  & age2014 > 61 & age2014 < 65 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 58 & age2014 < 62, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("59-61") append 
mean diedpanel [pweight=pwgt] if target2==1  & age2014 > 58 & age2014 < 62 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 55 & age2014 < 59, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("56-58") append 
mean diedpanel [pweight=pwgt] if target2==1  & age2014 > 55 & age2014 < 59 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 52 & age2014 < 56, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("53-55") append 
mean diedpanel [pweight=pwgt] if target2==1  & age2014 > 52 & age2014 < 56 & posttreated==1

*Broader bins:

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 18 & age2014 < 30, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("19-29") append 
mean diedpanel [pweight=pwgt] if target2==1 & age2014 > 18 & age2014 < 30 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 29 & age2014 < 40, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("30-39") append 
mean diedpanel [pweight=pwgt] if target2==1  & age2014 > 29 & age2014 < 40 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 39 & age2014 < 50, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("40-49") append 
mean diedpanel [pweight=pwgt] if target2==1 & age2014 > 39 & age2014 < 50 & posttreated==1

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target2==1 & age2014 > 49 & age2014 < 55, absorb(st dyear year) cluster(st)
outreg2 using $output/ages_3.xls, ctitle("50-54") append 
mean diedpanel [pweight=pwgt] if target2==1 & age2014 > 49 & age2014 < 55 & posttreated==1

*************************************************************************************************
***Placebo analyses
**Placebo--medicare eligible at time of the ACA
**Low income elderly
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if age2014 > 64 & target1b==1, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, append 

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if age2014 > 64 & target1b==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs.xls, append bdec(10)

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if age2014 > 64, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff_controls.xls, append bdec(10)

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if age2014 > 64, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs.xls, append bdec(10)

**Placebo---high income
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if povpi > 399, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, append 

cap noi reghdfe diedpanel posttreated sameyear $controls [pweight=pwgt] if povpi > 399, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff_controls.xls, append bdec(10)

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if povpi > 399 & age2014 > 54 & age2014 < 65, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs.xls, append bdec(10)

**Placebo---high income and age 55-64
cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if povpi > 399 & age2014 > 54 & age2014 < 65, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff.xls, append 

**Target sample + Age 55-64 in 2014 + not insured
cap noi reghdfe diedpanel posttreated sameyear $controls [pweight=pwgt] if target1==1 & age2014 > 54 & noinsure==1, absorb(st dyear year) cluster(st)
outreg2 using $output/diffdiff_controls.xls, append bdec(10)

*************************************************************************************
**Linear pre-trend Version 1 

gen timetoexp=0
replace timetoexp=dyear-2017 if treated17==1
replace timetoexp=dyear-2016 if treated16==1
replace timetoexp=dyear-2015 if treated15==1
replace timetoexp=dyear-2014 if treated14==1

cap noi reghdfe diedpanel event4 event3 event2 event1 timetoexp sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Pre-Trend 1") append

cap noi reghdfe diedpanel posttreated sameyear [pweight=pwgt] if target1==1 & age2014 > 54, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("Pre-Trend 1") append

*****************************************************************************************
***State pre-trends use 2-step procedure from AGB guide
preserve
keep if target1==1 & inrange(age2014,55,64)

qui: tab st, gen(stdum)
su stdum*
forval x=2/46{
  gen trend`x'=dyear-2008
  replace trend`x'=0 if stdum`x'==0
}
cap noi reghdfe diedpanel trend2-trend46 sameyear ///
[pweight=pwgt] if dyear<2014, absorb(st dyear year) 
predict yhat, xb
gen e=diedpanel-yhat

cap noi reghdfe e posttreated sameyear [pweight=pwgt], absorb(st dyear year) cluster(st)
outreg2 using $output/statepretrend.xls, ctitle("State pretrend") append 

cap noi reghdfe e event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt], absorb(st dyear year) cluster(st)
outreg2 using $output/statepretrend.xls, ctitle("State pretrend") append 
restore
*


************************************************************************
**Within treated states diff diff - 400%+ FPL as control group

gen fpl400=0
replace fpl400=1 if povpi > 399 & povpi!=.

forval n=1/4{
gen event400`n'=fpl400*event`n'
}

forval n=1/5{
gen event400m`n'=fpl400*eventm`n'
}

gen target1c=0
replace target1c=1 if povpi < 139
replace target1c=1 if povpi > 399 & povpi!=.
replace target1c=1 if schl < 16
replace target1c=0 if ssi > 0
replace target1c=0 if cit=="5"
replace target1c=0 if age2014 < 55
replace target1c=0 if age2014 > 64
replace target1c=0 if povpi > 999

gen posttreatedc=posttreated
replace posttreatedb=0 if fpl400==1

cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 event4c event3c event2c event1c eventm1c eventm2c eventm3c eventm4c eventm5c sameyear fpl400 fpl400##st fpl400##dyear [pweight=pwgt] if target1c==1 & everexpand==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("DDD-400%FPL As Controls") append 

cap noi reghdfe diedpanel posttreated posttreatedc sameyear fpl400 fpl##st fpl##dyear  [pweight=pwgt] if target1c==1 & age2014 > 54 & everexpand==1, absorb(st dyear year) cluster(st)
outreg2 using $output/panelregs_qje.xls, ctitle("DDD-400%FPL As Controls") append 

***********************************************************************
***Cluster by Census division
gen cendiv=.
replace cendiv=1 if inlist(st, 9,23,25,33,44,50)
replace cendiv=2 if inlist(st, 34,36,42)
replace cendiv=3 if inlist(st,18,17,26,39,55)
replace cendiv=4 if inlist(st,19,20,27,29,31,38,46)
replace cendiv=5 if inlist(st,10,11,12,13,24,37,45,51,54)
replace cendiv=6 if inlist(st,1,21,28,47)
replace cendiv=7 if inlist(st,5,22,40,48)
replace cendiv=8 if inlist(st,4,8,16,35,30,49,32,56)
replace cendiv=9 if inlist(st,2,6,15,41,53)

sysdir
adopath ++ .

set seed 1028455859
cap noi reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear [pweight=pwgt] if target1==1 & age2014 > 54 , absorb(st dyear year) cluster(cendiv)
outreg2 using $output/panelregs_qje.xls, ctitle("Cluster Census Division") append 
boottest event4=0, weight(webb)
boottest event3=0, weight(webb) 
boottest event2=0, weight(webb) 
boottest event1=0, weight(webb) 
boottest eventm1=0, weight(webb) 
boottest eventm2=0, weight(webb)
boottest eventm3=0, weight(webb) 
boottest eventm4=0, weight(webb) 
boottest eventm5=0, weight(webb) 

cap noi reghdfe diedpanel posttreated sameyear dyeardum2-dyeardum10 yeardum2-yeardum6 [pweight=pwgt], absorb(st) cluster(cendiv)
boottest posttreated=0, weight(webb)


***************************************************************************
*AGB DECOMPOSITION
* Based on (and adapted from) bacondecomp (Goodman-Bacon, Goldring, Nichols 15sep2019)

**AGB decomp
use $data/panelfinal, clear
keep if target1==1 & inrange(age2014,55,64)

*quick check that no missing vals on any reg vars
*su diedpanel posttreated sameyear dyear st year pwgt 

*DD model
gen g=1 if treated14==1
replace g=2 if treated15==1
replace g=3 if treated16==1
replace g=4 if treated17==1
replace g=5 if treated14==0 & treated15==0 & treated16==0 & treated17==0
lab def g 5 "Never", modify

tab year, gen(yr)
reg diedpanel posttreated sameyear yr2-yr6 i.st i.dyear [pweight=pwgt], cluster(st)
scalar DDest=_b[posttreated]
mat origb=e(b)
mat origv=e(V)

*Frisch-Waugh-Lovell regression
reg posttreated sameyear yr2-yr6 i.st i.dyear [pweight=pwgt], cluster(st)
predict double p if e(sample)

*collapse x's and p to the group/year level
loc xglist
foreach x in sameyear yr2 yr3 yr4 yr5 yr6 p {
  bys g dyear: g double sumg=sum(`x'*pwgt)
  by g dyear: g double sumwt=sum(pwgt/(`x'<.))
  by g dyear: g double g`x'=sumg[_N]/sumwt[_N]
  cap drop sumg sumwt
  if "`x'"!="p" loc xglist `xglist' g`x'
}
di "`xglist'"

*comparisons across all timing groups
foreach x in T C {
  g str1 `x'=""
}
foreach x in S B R2 cgroup {
  g byte `x'=.
}
la var T  "Bacon decomp 2x2 treated group"
la var C "Bacon decomp 2x2 control group"
la var S "Bacon decomp 2x2 weight"
la var B "Bacon decomp 2x2 coefficient"
la var R2 "Bacon decomp 2x2 R-squared"
cap la var g "Bacon decomp timing group"
la var cgroup "Bacon decomp aggregation group"

tab g posttreated, m

loc index 1
forv it=2/5{
  local jtmax=`it'-1
  forv jt=1/`jtmax'{
    di "Loop it=`it', jt=`jt'"
    loc itstring="`:lab g `it''"
    loc jtstring="`:lab g `jt''"
    g byte samp=((g==`it')|(g==`jt'))
    su samp [aw=pwgt], mean
    scalar share=r(mean)

    *VD so x's not partialled out
    *get dyad variance
    reg posttreated i.st i.dyear [pweight=pwgt] if samp, cluster(st)
    predict double dum if e(sample)
    gen Dtilde=posttreated-dum if e(sample)
    drop dum
    sum Dtilde [aw=pwgt]
    scalar VD=((r(N)-1)/r(N))*r(Var)


    *partials FE out of the group-level x's, not individ-level hence gvar
    local XXlist
    foreach var of varlist sameyear yr2-yr6 {
	reg g`var' i.st i.dyear [pweight=pwgt] if samp, cluster(st)
	predict double dum if e(sample)
	gen XX`var'=g`var'-dum if e(sample)
	drop dum
	local XXlist `XXlist' XX`var'
     }
    cap drop pgjtilde
    reg Dtilde `XXlist' [aw=pwgt]
    scalar Rsq=e(r2)
    predict double pgjtilde if e(sample)
    drop `XXlist'
    *partialled out p in the dyad
    cap drop ptilde
    reg gp i.st i.dyear [pweight=pwgt] if samp, cluster(st)
    predict double dum if e(sample)
    gen ptilde=gp-dum if e(sample)
    drop dum
    gen double dp=pgjtilde-ptilde if samp
    *get variance of pg tilde in the dyad
    sum dp [aw=pwgt]
    scalar Vdp=((r(N)-1)/r(N))*r(Var)
    *weight: this is the variance of dtilde in the dyad
    scalar finals=(scalar(share))^2*((1-scalar(Rsq))*scalar(VD) + scalar(Vdp))

    *get the proper X-adjusted dyad coef
    reg diedpanel posttreated `xglist' i.st i.dyear [pweight=pwgt] if samp, cluster(st)
    scalar BD=_b[posttreated]
    reg diedpanel dp i.st i.dyear [pweight=pwgt] if samp, cluster(st)
    scalar Bb=_b[dp]
    *The dyad "Beta" combines the proper controlled one and a term for how wrong the "FWL" coef is
    scalar Beta=((1-scalar(Rsq))*scalar(VD)*scalar(BD) + ///
    scalar(Vdp)*scalar(Bb))/((1-scalar(Rsq))*scalar(VD) + scalar(Vdp))
  

    qui {
      replace T="`itstring'" in `index'
      replace C="`jtstring'" in `index'
      replace S=scalar(finals) in `index'
      replace B=scalar(Beta) in `index'
      replace R2=scalar(Rsq) in `index'
      replace cgroup=1*(!inlist("`itstring'","Never")&!inlist("`jtstring'","Never"))+ ///
      3*inlist("Never","`itstring'","`jtstring'") in `index'
      loc index=`index'+1
    }
    foreach name in samp Dtilde dp {
      cap drop `name'
    }
  }
}
lab def cgroup 1 "Timing groups" 2 "Always treated vs timing" 3 "Never treated vs timing" ///
4 "Within" 
la val cgroup cgroup

*within part
reg diedpanel i.g##i.dyear p i.st [pweight=pwgt], cluster(st)
scalar Bw=-_b[p]
reg p i.g##i.dyear i.st [pweight=pwgt], cluster(st)
predict double dum if e(sample)
gen pw=p-dum if e(sample)
drop dum
sum pw [aw=pwgt]
scalar Vw=((r(N)-1)/r(N))*r(Var)
replace T="Within" in `index'
replace C="" in `index'
replace S=scalar(Vw) in `index'
replace B=scalar(Bw) in `index'
replace cgroup=4 in `index'

*rescale weights
su S if T~="Within", mean
scalar tb=r(sum)
gen double sg=S/scalar(tb) if T!="Within"
su S, mean
scalar totals=r(sum)
g double omega=1-scalar(tb)/scalar(totals)
su omega, mean
scalar O=r(sum)
gen double DD=sg*(1-scalar(O))*B*(T~="Within")+scalar(O)*B*(T=="Within") if !mi(S)
gen double wgt=S/scalar(totals) if !mi(S)

*post estimates
forv ix=1/`index'{
  mat sumcalc=nullmat(sumcalc)\(cgroup[`ix'],B[`ix'],wgt[`ix'])
  if cgroup[`ix']==4 loc winest=B[`ix']
  if cgroup[`ix']==4 loc winwgt=wgt[`ix']
  if cgroup[`ix']==4 loc cvcwgt=wgt[`ix']
  mat postb=nullmat(postb), B[`ix']
  mat postv=nullmat(postv), wgt[`ix']
  mat output1=nullmat(output1), B[`ix']
  mat output2=nullmat(output2), wgt[`ix']
  loc rname `rname' "`=T[`ix']'`=cond(C[`ix']=="","","_")'`=C[`ix']'"
}
mat li sumcalc
mata:st_matrix("summary2",panelsum(st_matrix("sumcalc")[.,3],panelsetup(st_matrix("sumcalc"),1)))
mata:st_matrix("summary1",panelsum(st_matrix("sumcalc")[.,2],st_matrix("sumcalc")[.,3],panelsetup(st_matrix("sumcalc"),1)):/panelsum(st_matrix("sumcalc")[.,3],panelsetup(st_matrix("sumcalc"),1)))
mata:st_matrix("summary3",panelsum(st_matrix("sumcalc")[.,1],panelsetup(st_matrix("sumcalc"),1)):/panelsum(st_matrix("sumcalc")[.,1]:<.,panelsetup(st_matrix("sumcalc"),1)))
mat summary=summary1, summary2
*row 1=timing groups, row 2=never v timing, row 3=within
mat li summary

*component estimates and weights
mat li postb
mat li postv

**************************************************************************************
*Sun & Abraham estimation, based on (and adapted from) replication code for Sun & Abraham (2020)

use $data/panelfinal, clear
keep if target1==1 & inrange(age2014,55,64)

global event eventm5 eventm4 eventm3 eventm2 eventm1 event1 event2 event3 event4

**two-way FE
reghdfe diedpanel $event sameyear [pweight=pwgt], absorb(st dyear year) cluster(st)

**S&A interaction-weighted estimator
svyset [pweight=pwgt]
*rename and reorder event time, create timing group dummies
rename eventm5 evttime1
rename eventm4 evttime2
rename eventm3 evttime3
rename eventm2 evttime4
rename eventm1 evttime5
rename event0 evttime6
rename event1 evttime7
rename event2 evttime8
rename event3 evttime9
rename event4 evttime10
order evttime1 evttime2 evttime3 evttime4 evttime5 evttime6 evttime7 evttime8 ///
evttime9 evttime10
gen expyr1=treated14==1
gen expyr2=treated15==1
gen expyr3=treated16==1
gen expyr4=treated17==1

*get variance from estimating share of cohorts
preserve
foreach vv of varlist evttime* {
  replace `vv'=. if treated14==0 & treated15==0 & treated16==0 & treated17==0 
}
estimates clear
forval y=1/4{
  svy: reg expyr`y' evttime1-evttime9, nocons
  estimates store regexp`y'
  mat z`y'=e(b)
}
suest regexp1 regexp2 regexp3 regexp4
mat V=e(V)
mat Vcov=V
mat list Vcov
restore

*saturated specification
foreach vv of varlist evttime* {
  forval i=1(1)4{
    gen `vv'_`i'=`vv'*(expyr`i'==1)
  }
}
global rhs_rel_year_i ""
global rhs_rel_year_i "$rhs_rel_year_i evttime1_1 evttime2_1 evttime3_1 evttime4_1 evttime5_1 evttime7_1 evttime8_1 evttime9_1 evttime10_1"
global rhs_rel_year_i "$rhs_rel_year_i evttime1_2 evttime2_2 evttime3_2 evttime4_2 evttime5_2 evttime7_2 evttime8_2 evttime9_2"
global rhs_rel_year_i "$rhs_rel_year_i evttime1_3 evttime2_3 evttime3_3 evttime4_3 evttime5_3 evttime7_3 evttime8_3"
global rhs_rel_year_i "$rhs_rel_year_i evttime1_4 evttime2_4 evttime3_4 evttime4_4 evttime5_4 evttime7_4"

reghdfe diedpanel $rhs_rel_year_i sameyear [pweight=pwgt], absorb(st dyear year) cluster(st)
mat b=e(b)
mat list b
mat V=vecdiag(e(V))
mat b_1=(b[.,1..5],0,b[.,6..9])
mat V_1=(V[.,1..5],0,V[.,6..9])
mat list b_1
mat list V_1
mat b_2=(b[.,10..14],0,b[.,15..17],.)
mat V_2=(V[.,10..14],0,V[.,15..17],.)
mat list b_2
mat list V_2
mat b_3=(b[.,18..22],0,b[.,23..24],.,.)
mat V_3=(V[.,18..22],0,V[.,23..24],.,.)
mat list b_3
mat list V_3
mat b_4=(b[.,25..29],0,b[.,30..30],.,.,.)
mat V_4=(V[.,25..29],0,V[.,30..30],.,.,.)
mat list b_4
mat list V_4

*take weighted average of estimates from the saturated specification to form IW estimates
mat biw_long=b

*evttime1
matselrc biw_long tempb, c(1,10,18,25)
matselrc Vcov tempVcov, r(1,10,19,28) c(1,10,19,28)
mat list tempb
mat list tempVcov
mat temp=tempb*tempVcov*tempb'
lincom z1[1,1]*evttime1_1+z2[1,1]*evttime1_2+z3[1,1]*evttime1_3+z4[1,1]*evttime1_4
mat b=r(estimate)
mat V=(r(se)^2+temp)

*evttime2
matselrc biw_long tempb, c(2,11,19,26)
matselrc Vcov tempVcov, r(2,11,20,29) c(2,11,20,29)
mat list tempb
mat list tempVcov
mat temp=tempb*tempVcov*tempb'
lincom z1[1,2]*evttime2_1+z2[1,2]*evttime2_2+z3[1,2]*evttime2_3+z4[1,2]*evttime2_4
mat b=b\r(estimate)
mat V=V\(r(se)^2+temp)

*evttime3
matselrc biw_long tempb, c(3,12,20,27)
matselrc Vcov tempVcov, r(3,12,21,30) c(3,12,21,30)
mat list tempb
mat list tempVcov
mat temp=tempb*tempVcov*tempb'
lincom z1[1,3]*evttime3_1+z2[1,3]*evttime3_2+z3[1,3]*evttime3_3+z4[1,3]*evttime3_4
mat b=b\r(estimate)
mat V=V\(r(se)^2+temp)

*evttime4
matselrc biw_long tempb, c(4,13,21,28)
matselrc Vcov tempVcov, r(4,13,22,31) c(4,13,22,31)
mat list tempb
mat list tempVcov
mat temp=tempb*tempVcov*tempb'
lincom z1[1,4]*evttime4_1+z2[1,4]*evttime4_2+z3[1,4]*evttime4_3+z4[1,4]*evttime4_4
mat b=b\r(estimate)
mat V=V\(r(se)^2+temp)

*evttime5
matselrc biw_long tempb, c(5,14,22,29)
matselrc Vcov tempVcov, r(5,14,23,32) c(5,14,23,32)
mat list tempb
mat list tempVcov
mat temp=tempb*tempVcov*tempb'
lincom z1[1,5]*evttime5_1+z2[1,5]*evttime5_2+z3[1,5]*evttime5_3+z4[1,5]*evttime5_4
mat b=b\r(estimate)
mat V=V\(r(se)^2+temp)

*evttime6
mat b=b\0
mat V=V\0

*evttime7
matselrc biw_long tempb, c(6,15,23,30)
matselrc Vcov tempVcov, r(7,16,25,34) c(7,16,25,34)
mat list tempb
mat list tempVcov
mat temp=tempb*tempVcov*tempb'
lincom z1[1,7]*evttime7_1+z2[1,7]*evttime7_2+z3[1,7]*evttime7_3+z4[1,7]*evttime7_4
mat b=b\r(estimate)
mat V=V\(r(se)^2+temp)

*evttime8
matselrc biw_long tempb, c(7,16,24)
matselrc Vcov tempVcov, r(8,17,26) c(8,17,26)
mat list tempb
mat list tempVcov
mat temp=tempb*tempVcov*tempb'
lincom z1[1,8]*evttime8_1+z2[1,8]*evttime8_2+z3[1,8]*evttime8_3
mat b=b\r(estimate)
mat V=V\(r(se)^2+temp)

*evttime9
matselrc biw_long tempb, c(8,17)
matselrc Vcov tempVcov, r(9,18) c(9,18)
mat list tempb
mat list tempVcov
mat temp=tempb*tempVcov*tempb'
lincom z1[1,9]*evttime9_1+z2[1,9]*evttime9_2
mat b=b\r(estimate)
mat V=V\(r(se)^2+temp)

*evttime10
lincom evttime10_1
mat b=b\r(estimate)
mat V=V\r(se)^2

mat biw=b,V,b_1',V_1',b_2',V_2',b_3',V_3',b_4',V_4'
mat list biw

clear 
global n_obs=10
set obs $n_obs

gen y=_n-6
svmat biw

gen sd2=sqrt(biw2)

list y biw1 sd2
export excel using $output/sunabraham.xlsx, replace

****************************************************************************************************
*Conley (2008) spatial standard errors
use $data/panelfinal, clear
keep if target1==1 & inrange(age2014,55,64)
keep diedpanel posttreated sameyear st year dyear centroidlat centroidlon pwgt countyfips event*
su year dyear centroidlat centroidlon countyfips 
gen ob=1
tab year, gen(yr)
collapse diedpanel posttreated sameyear yr2 yr3 yr4 yr5 yr6 event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 ///
eventm5 (sum) pop=ob [pweight=pwgt], by(countyfips centroidlat centroidlon dyear st)
merge m:1 countyfips using $data/fixcentroid, update 

reghdfe diedpanel posttreated sameyear yr2 yr3 yr4 yr5 yr6 [aweight=pop], absorb(st dyear) cluster(st)

acreg diedpanel posttreated sameyear yr2 yr3 yr4 yr5 yr6 [fweight=pop], id(st) time(dyear) spatial latitude(centroidlat) longitude(centroidlon) ///
pfe1(st) pfe2(dyear) lag(10) dist(500) hac bartlett

reghdfe diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear yr2 yr3 yr4 yr5 yr6 [fweight=pop], ///
absorb(st dyear) cluster(st)

acreg diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear yr2 yr3 yr4 yr5 yr6 [fweight=pop], ///
id(st) time(dyear) spatial latitude(centroidlat) longitude(centroidlon) pfe1(st) pfe2(dyear) lag(10) dist(500) hac bartlett

*****************************************************************************************************
*Nonlinear models
use $data/panelfinalelig, clear
keep if target1==1 & inrange(age2014,55,64)
save $data/nonlinear, replace

*Logit
logit diedpanel posttreated sameyear i.st i.dyear i.year [pweight=pwgt], cluster(st)
outreg2 posttreated using $output/logit.xls, replace
margins, dydx(posttreated)

logit diedpanel event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 sameyear i.st i.year i.dyear [pweight=pwgt], cluster(st)
outreg2 event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 using $output/logit.xls, append
margins, dydx(event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5)

*Cox
use $data/nonlinear, clear
keep dyear studyid deathyear died year dyear diedpanel sameyear st pwgt posttreated event* ///
treated*

egen idvar=concat(year studyid)
stset dyear [pweight=pwgt], failure(diedpanel) id(idvar) origin(time year)

qui: tab st, gen(sta)
qui: tab dyear, gen(dyr)
qui: tab year, gen(yr)


stcox posttreated yr1-yr5 sta1-sta45 sameyear dyr1-dyr9, vce(cluster st)
outreg2 posttreated using $output/cox.xls, replace

stcox event4 event3 event2 event1 eventm1 eventm2 eventm3 eventm4 eventm5 yr1-yr5 sta1-sta45 sameyear dyr1-dyr9, vce(cluster st)
outreg2 posttreated using $output/cox.xls, append

******************************************************************************************************
*Placebo simulations for randomly assigned expansion statusu during pre-ACA years 2004-2013 
forval y=2004/2013{
	use $data/pers`y', clear
	keep age pwgt st povpi died deathyear schl ssi cit
	gen year=`y'
	destring st, replace
	drop if inlist(st,11,10,25,36,50)
	gen age2010=2010-year+age
	keep if inrange(age2010,55,64)
	destring schl, replace
	gen target1=0
	replace target1=1 if povpi<139
	if `y'<2008{
	  replace target1=1 if schl<9
	}
	if `y'>=2008{
	  replace target1=1 if schl < 16
	}
	replace target1=0 if ssi > 0
	replace target1=0 if cit=="5"
	keep if target1==1
	di "Year `y'"   
	count
	save $data/sim`y', replace
}

*correct for 2004 sample being approx 30% of the size of other years
*expand by 3
use $data/sim2004, clear
expand 3
replace pwgt=pwgt/3
save $data/sim2004, replace

forval t=2004/2009 {
  use $data/sim`t', clear
  gen obnum=_n
  egen id=concat(year obnum)
  save $data/sim`t', replace

  forval y=`t'/2013{
    use $data/sim`t', clear
    gen dyear=`y'
    gen diedpanel=0
    replace diedpanel=1 if deathyear==`y'
    replace diedpanel=. if deathyear < `y'
    gen sameyear=0
    replace sameyear=1 if dyear==year
    save $data/panel`t'`y', replace
  }

  use $data/panel`t'`t', clear
  forval y=`t'/2013{
    if `y'!=`t'{
      append using $data/panel`t'`y'
    }
  }
  save $data/panel`t'all, replace
}

use $data/panel2004all, clear
forval y=2005/2009{
  append using $data/panel`y'all
}
save $data/panelsim, replace
su
bysort dyear: su

global numexpand2010=21
global numexpand2011=3
global numexpand2012=2
global numexpand2013=1

putexcel set "$output/placebosim.xlsx", replace
putexcel A1=("coef")
putexcel B1=("se")
putexcel C1=("df")

timer on 1
local max_num=10000
local max_row=`max_num'
local row=2
forval y=1/`max_num'{
  di "`y' iteration"

  use $data/panelsim, clear
  local rand_seed=81635+`y'
  set seed `rand_seed'

  bysort st: gen dum=runiform() if _n==1
  *su dum
  bysort st: egen randvar=max(dum)
  *su randvar
  egen rank=group(randvar)
  *su rank
  gen posttreated=0
  replace posttreated=1 if rank<=$numexpand2010 & dyear>=2010
  replace posttreated=1 if rank>$numexpand2010 & rank<=$numexpand2010+ ///
  $numexpand2011 & dyear>=2011
  replace posttreated=1 if rank>$numexpand2010+$numexpand2011 & rank<=$numexpand2010+ ///
  $numexpand2011+$numexpand2012 & dyear>=2012
  replace posttreated=1 if rank>$numexpand2010+$numexpand2011+$numexpand2012 & ///
  rank<=$numexpand2010+$numexpand2011+$numexpand2012+$numexpand2013 & dyear>=2013
  *bysort dyear: su posttreated


  reghdfe diedpanel posttreated  [pweight=pwgt], absorb(st dyear year sameyear) cluster(st)
  putexcel A`row'=(_b[posttreated])
  putexcel B`row'=(_se[posttreated])
  putexcel C`row'=(`e(df_r)')

  di "PERCENT COMPLETE"
  di ((`row'-1)/`max_row')*100
  local row=`row'+1
}
timer off 1 
timer list
*/

***addition on 2/28/20 to graph output
import excel "$output/placebosim.xlsx", sheet("Sheet1") firstrow clear
describe
list in 1/10

su coef, d
local ci_low=r(p5)
local ci_high=r(p95)
sort coef

graph twoway (histogram coef, fcolor(gs12) lcolor(black)) (scatteri 0 -0.00132 525 -0.00132, c(l) m(i) lcolor(red) lpattern(dash)) ///
(scatteri 0 `ci_low' 525 `ci_low', c(l) m(i) lcolor(blue)) ///
(scatteri 0 `ci_high' 525 `ci_high', c(l) m(i) lcolor(blue)), ///
legend(off) graphregion(color(white)) 
graph export "$output/placebo_coef.jpg", replace as(jpg)

    
gen tstat=coef/se

su tstat, d
local ci_low=r(p5)
local ci_high=r(p95)

graph twoway (histogram tstat, fcolor(gs12) lcolor(black)) (scatteri 0 -2.656 .325 -2.656, c(l) m(i) lcolor(red) lpattern(dash)) ///
(scatteri 0 `ci_low' .325 `ci_low', c(l) m(i) lcolor(blue)) ///
(scatteri 0 `ci_high' .325 `ci_high', c(l) m(i) lcolor(blue)), legend(off) graphregion(color(white))
graph export "$output/placebo_tstat.jpg", replace as(jpg)







    
